source("~/Documents/GitHub/VivID_Epi/analyses/00-functions.R")
library(tidyverse)
library(sf)
library(srvyr) #wrap the survey package in dplyr syntax
#......................
# Import Data
#......................
load("~/Documents/GitHub/VivID_Epi/data/vividepi_raw.rda")
dt <- merge_pr_plsmdm_gemtdt(pr = arpr, plsmdm = panplasmpcrres, ge = ge)Steve has asked me to produce maps and basic correlation plots of Pf versus Pv
prev_summarizer <- function(data, maplvl, plsmdmspec, sfobj){
# catch error
# if( !( any(colnames(x) %in% c("hv001")) & any(colnames(x) %in% c("hv005")) ) ){
# stop("Must include cluster and cluster weight. need to check if you are using pr or hiv sample weights")
# }
# rlang
maplvl <- enquo(maplvl)
plsmdmspec <- enquo(plsmdmspec)
# clusters are weighted (each individual has same weight in cluster)
ret <- data %>%
srvyr::as_survey_design(ids = hv001, weights = hv005) %>%
dplyr::group_by(!!maplvl) %>%
dplyr::summarise(plsmdn = srvyr::survey_total(hv001),
plsmd = srvyr::survey_mean(!!plsmdmspec, na.rm = T, vartype = c("se", "ci"), level = 0.95)) %>%
dplyr::left_join(., sfobj) # attach spatial data, let R figure out the common var
# return
return(ret)
}
pfldhprov <- prev_summarizer(data = dt, maplvl = adm1name, plsmdmspec = pfldh, sfobj = DRCprov) %>%
dplyr::mutate(plsmdmspec = "pfldh", maplvl = "adm1name")
pv18sprov <- prev_summarizer(data = dt, maplvl = adm1name, plsmdmspec = pv18s, sfobj = DRCprov) %>%
dplyr::mutate(plsmdmspec = "pv18s", maplvl = "adm1name")
pfldhclust <- prev_summarizer(data = dt, maplvl = hv001, plsmdmspec = pfldh, sfobj = ge) %>%
dplyr::mutate(plsmdmspec = "pfldh", maplvl = "hv001")
pv18sclust <- prev_summarizer(data = dt, maplvl = hv001, plsmdmspec = pv18s, sfobj = ge) %>%
dplyr::mutate(plsmdmspec = "pv18s", maplvl = "hv001")
mp <- dplyr::bind_rows(pfldhprov, pv18sprov, pfldhclust, pv18sclust) %>%
dplyr::group_by(plsmdmspec, maplvl) %>%
tidyr::nest()
mp$data <- sapply(list(pfldhprov, pv18sprov, pfldhclust, pv18sclust), function(x) return(x))
mapplotter <- function(data, maplvl, plsmdmspec){
ret <- ggplot() +
geom_sf(data = DRCprov) +
ggtitle(paste(plsmdmspec)) +
theme(plot.title = element_text(famil = "Arial", face = "bold", hjust = 0.5, size = 13),
legend.position = "bottom",
legend.title = element_text(famil = "Arial", face = "bold", vjust = 0.85, size = 12),
legend.text = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, angle = 90, size = 10),
axis.text = element_blank(),
axis.line = element_blank(),
panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent"),
panel.grid = element_blank(),
panel.border = element_blank())
if(maplvl == "adm1name"){
ret <- ret + geom_sf(data = data, aes(fill = plsmd)) +
scale_fill_gradient2("Prevalence", low = "#0000FF", mid = "#FFEC00", high = "#FF0000") +
coord_sf(datum=NA) # to get rid of gridlines
} else if(maplvl == "hv001"){
ret <- ret + geom_sf(data = data, aes(fill = plsmd, colour = plsmd, size = plsmdn), alpha = 0.8) +
scale_color_gradient2("Prevalence", low = "#0000FF", mid = "#FFEC00", high = "#FF0000") +
scale_size(guide = 'none') + scale_fill_continuous(guide = 'none') +
coord_sf(datum=NA) # to get rid of gridlines
} else {
stop("maplvl is not in the options for this function")
}
return(ret)
}
prevmaps <- pmap(mp, mapplotter)gridExtra::grid.arrange(prevmaps[[1]],
prevmaps[[3]],
ncol=2, top=grid::textGrob("Pf Prevalence by Province and Cluster in CD2013 DHS", gp=grid::gpar(fontsize=15, fontfamily = "Arial", fontface = "bold"))) gridExtra::grid.arrange(prevmaps[[2]],
prevmaps[[4]],
ncol=2, top=grid::textGrob("Pf Prevalence by Province and Cluster in CD2013 DHS", gp=grid::gpar(fontsize=15, fontfamily = "Arial", fontface = "bold"))) gridExtra::grid.arrange(prevmaps[[1]],
prevmaps[[2]],
prevmaps[[3]],
prevmaps[[4]],
ncol=2, nrow=2, top=grid::textGrob("Prevalence by Species and Province/Cluster in CD2013 DHS", gp=grid::gpar(fontsize=15, fontfamily = "Arial", fontface = "bold"))) Points are scaled by size of province (and color of province). Loess curve in red.
pfprov <- mp$data[[1]]
colnames(pfprov)[4] <- "Pf_prev"
pvprov <- mp$data[[2]]
colnames(pvprov)[4] <- "Pv_prev"
plotObj <- dplyr::bind_cols(pfprov, pvprov) %>%
ggplot(aes(x=Pf_prev, y=Pv_prev, size = plsmdn)) +
geom_point(aes(colour = adm1name)) +
geom_smooth(method="loess", se=F, colour = "red") +
ggtitle("Pfalciparum versus Pvivax Province Prevalence") + ylab("Pv Prevalence") + xlab("Pf Prevalence") +
theme(plot.title = element_text(famil = "Arial", face = "bold", hjust = 0.5, size = 18),
legend.position = "none",
axis.text.x = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, angle = 90, size = 10),
axis.text.y = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, size = 10),
axis.ticks = element_blank(),
axis.line = element_line(colour = "black"),
panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent"),
panel.grid = element_blank(),
panel.border = element_blank())
plotly::ggplotly(plotObj)Points are scaled by size of province (and color of province). Loess curve in red. Cluster filtered for both 0s.
pfprov <- mp$data[[3]]
colnames(pfprov)[4] <- "Pf_prev"
pvprov <- mp$data[[4]]
colnames(pvprov)[4] <- "Pv_prev"
plotObj <- dplyr::bind_cols(pfprov, pvprov) %>%
dplyr::filter(!(Pf_prev == 0 & Pv_prev == 0)) %>%
ggplot(aes(x=Pf_prev, y=Pv_prev, size = plsmdn)) +
geom_point(aes(colour = hv001)) +
geom_smooth(method="loess", se=F, colour = "red") +
ggtitle("Pfalciparum versus Pvivax Province Prevalence") + ylab("Pv Prevalence") + xlab("Pf Prevalence") +
theme(plot.title = element_text(famil = "Arial", face = "bold", hjust = 0.5, size = 18),
legend.position = "none",
axis.text.x = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, angle = 90, size = 10),
axis.text.y = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, size = 10),
axis.ticks = element_blank(),
axis.line = element_line(colour = "black"),
panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent"),
panel.grid = element_blank(),
panel.border = element_blank())
plotly::ggplotly(plotObj)